/* Finite element subroutines: mesh data structures
   Jan Stebel
*/

#ifndef _MESH_H__
#define _MESH_H__

#include <fstream>
#include <string>
#include <iostream>
#include <map>
#include <vector>
#include <set>

#include "ad.h"
#include "fe_types.h"

#define GMSH_LINE             1
#define GMSH_TRIANGLE         2
#define GMSH_TETRAHEDRON      4
#define GMSH_QUAD_LINE        8
#define GMSH_QUAD_TRIANGLE    9
#define GMSH_QUAD_TETRAHEDRON 11
#define GMSH_POINT            15

#define VTK_LINE               3
#define VTK_TRIANGLE           5
#define VTK_TETRA              10
#define VTK_QUAD_EDGE     21
#define VTK_QUAD_TRIANGLE 22
#define VTK_QUAD_TETRA    24

/** VTK ordering of nodes differs from GMSH */
const int VTK_QUAD_TETRA_NODE_ORDER[] = { 0, 1, 2, 3, 4, 5, 6, 7, 9, 8 };


using namespace std;

typedef vector<int> ChildParentListType;

/** Structure for elementary geometric entity */
class CElement
{
// protected:
//  inline void allocateNodes(int n) { nodes = new int[n]; nnodes = n; };
  
 public:
  int nphysical;  /**< Number of physical entity */
  int ngeom;      /**< Number of elementary geometrical entity */
  int npart;      /**< Number of mesh partition */
  int type;       /**< Element type */
  ChildParentListType children; /**< Child elements (e.g. faces of a volume) */
  ChildParentListType parents; /**< Parent elements (e.g. volumes adjacent to a face) */
  vector<int> nodes;    /**< Array of node numbers. */
  map<int, map<int, int> > parentNodes; /**< Node ordering of parent element */
  
  CElement() {};
  CElement(const CElement &e);
  ~CElement() { nodes.clear(); };
};


/** Abstract class for 3D mesh. */
class Mesh
{
 protected:
  vector<CElement> volumes, faces, edges, points;
  int nnodes, nlnodes;
  vector<int> lnmap;
  vector<int> node_partition;    /**< Partition numbers of mesh nodes. */
  map<int,set<int> > partition_nodes; /**< Node numbers belonging to partitions. */
  set<int> shared_nodes; /**< Nodes belonging to more than one partition. */
 
  vector<CVar>    coords;        /**< Node coordinates.                   */
  map<int,string> physicalNames; /**< Names of physical entities.         */
  map<string,int> physNameIds;   /**< Reverse array of physical entities. */
  
  int type;
  string fname;

  int addMidNode(int size, ...);
  
  
 public:
  /** Implicit constructor. */
  Mesh() {};
  
  /** Destructor. */
  virtual ~Mesh();
  
  /** Delete allocated arrays and structures. */
  virtual void destroy();
  
  /**
    Virtual method for reading mesh from file.
    @param fname File name.
    @param partition Read given part of the mesh (value 0 means: read the whole mesh).
  */
  void readGMSH (string ifname, int partition = 0);
  
  /**
    Virtual method for writing mesh to file.
    @param fname File name.
  */
  virtual void writeGMSH(string fname) = 0;
  
  int getType() { return type; };

  string fileName() { return fname; };

  void splitVolume(int v);
  
  void splitFace(int f);
  
  /**
    Deform the mesh coordinates according to given function.
    @param func Function performing deformation of comp-th coordinate
                of point located ad (x,y,z).
  */
  void         deform(CVar& func(double x, double y, double z, int comp));
  
  
  void splitBoundaryVolumes();
  
  
  virtual void findChildrenParents(bool create_new = false) = 0;
  signed int getFaceNparents   (int e) { return faces[e].parents.size(); };
  signed int getEdgeNparents   (int e) { return edges[e].parents.size(); };
  signed int getVolumeNchildren(int e) { return volumes[e].children.size(); };
  signed int getFaceNchildren  (int e) { return faces[e].children.size(); };
  int getFaceParentNo (int e, int n) { return faces[e].parents[n]; };
  int getEdgeParentNo (int e, int n) { return edges[e].parents[n]; };
  int getVolumeChildNo(int e, int n) { return volumes[e].children[n]; };
  int getFaceChildNo  (int e, int n) { return faces[e].children[n]; };
  int getFaceParentNodeNo(int e, int p, int n) { return faces[e].parentNodes[n][p]; };
  int getEdgeParentNodeNo(int e, int p, int n) { return edges[e].parentNodes[n][p]; };
  
  /**
    Return coordinate of mesh node.
    @param node  Node number.
    @param coord Coordinate number (0=x, 1=y, 2=z).
  */
  CVar& getCoord(int node, int coord);
  
  /**
    Set mesh node coordinate.
    @param node  Node number.
    @param coord Coordinate number.
    @param x     New coordinate value.
  */
  void  setCoord(int node, int coord, CVar x);
  
  /** Return number of physical entities. */
  int    getNphysicalNames()    { return physicalNames.size(); }
  
  /**
    Return physical name.
    @param n Number of physical entity.
  */
  string getPhysicalName(int n) { return physicalNames[n];     }
  
  /**
    Return number of physical entity given by name.
    @param name Name of physical entity.
  */
  int    getPhysNameId(string name) { return physNameIds[name]; }
  
  /** Iterator to beginning of physical names. */
  map<int,string>::iterator beginPhysicalNames() { return physicalNames.begin(); }
  
  /** Iterator to end of physical names. */
  map<int,string>::iterator endPhysicalNames() { return physicalNames.end(); }
  
  /**
    Calculate normal and tangent vectors on face.
    @param fid Face number.
    @param n   Normal vector.
    @param t1  First tangent vector.
    @param t2  Second tangent vector.
  */
  void   calculateFaceNormal(int fid, CVec &n, CVec &t1, CVec &t2);
  
  /**
    Calculate tangent vector on edge.
    @param eid Edge number.
    @param t   Tangent vector.
  */
  void   calculateEdgeTangent(int eid, CVec &t);
  
  /** Return number of nodes. */
  int getNnodes()     { return nnodes; }
  
  /**
   * Return number of nodes in a given partition.
   * @param part Partition number.
   */
  int getNnodes(int part)     { return partition_nodes[part].size(); }

  /** Return number of linear nodes. */
  int getNlnodes()    { return nlnodes; }
  
  /** Translate node no. within all nodes into linear node no. */
  int getNoOfLnode(int node) { return lnmap[node]; }
  
  /** Return number of point elements. */
  int getNpoints()     { return points.size(); }
  
  /** Return number of edges. */
  int getNedges()     { return edges.size(); }
  
  /** Return number of faces. */
  int getNfaces()     { return faces.size(); }
  
  /** Return number of volumes. */
  int getNvolumes()   { return volumes.size(); }

  /**
    Return physical entity of given volume.
    @param v Volume number.
  */
  int getVolumePhysicalNo(int v) { return volumes[v].nphysical; };
  
  /**
    Return physical entity of given face.
    @param f Face number.
  */
  int getFacePhysicalNo  (int f) { return faces[f].nphysical; };
  
  /**
    Return physical entity of given edge.
    @param e Edge number.
  */
  int getEdgePhysicalNo  (int e) { return edges[e].nphysical; };
  
  /**
    Return physical entity of given point.
    @param e Point number.
  */
  int getPointPhysicalNo  (int e) { return points[e].nphysical; };
  
  /**
    Return geometric entity of given volume.
    @param v Volume number.
  */
  int getVolumeGeomNo(int v) { return volumes[v].ngeom; };
  
  /**
    Return geometric entity of given face.
    @param f Face number.
  */
  int getFaceGeomNo  (int f) { return faces[f].ngeom; };
  
  /**
    Return geometric entity of given edge.
    @param e Edge number.
  */
  int getEdgeGeomNo  (int e) { return edges[e].ngeom; };
  
  /**
    Return geometric entity of given point.
    @param e Point number.
  */
  int getPointGeomNo  (int e) { return points[e].ngeom; };
  
  /**
    Return partition number of given volume.
    @param v Volume number.
  */
  int getVolumePartNo(int v) { return volumes[v].npart; };
  
  /**
    Return partition number of given face.
    @param f Face number.
  */
  int getFacePartNo  (int f) { return faces[f].npart; };
  
  /**
    Return partition number of given edge.
    @param e Edge number.
  */
  int getEdgePartNo  (int e) { return edges[e].npart; };
  
  /**
    Return partition number of given point.
    @param e Point number.
  */
  int getPointPartNo  (int e) { return points[e].npart; };
  
  /**
    Return partition number of given node.
    @param e Node number.
  */
  int getNodePartNo  (int e) { return node_partition[e]; };

  /**
    Return global index of given volume node.
    @param volume Volume number.
    @param node   Local node number.
  */
  int getVolumeNode(int volume, int node) { return volumes[volume].nodes[node]; };
  
  /**
    Return global index of given face node.
    @param face Face number.
    @param node Local node number.
  */
  int getFaceNode  (int face, int node) { return faces[face].nodes[node]; };
  
  /**
    Return global index of given edge node.
    @param edge Edge number.
    @param node Local node number.
  */
  int getEdgeNode  (int edge, int node) { return edges[edge].nodes[node]; };
  
  /**
    Return global index of point node.
    @param edge Point number.
  */
  int getPointNode  (int point) { return points[point].nodes[0]; };
  
  /**
   * Return set of nodes belonging to a partition.
   * @param part Partition number.
   */
  const set<int> &getPartitionNodes(int part) { return partition_nodes[part]; };

  /**
   * Return the set of shared nodes.
   */
  const set<int> &getSharedNodes() { return shared_nodes; };

  /** Return number of d.o.f. per volume. */
  virtual int getNvolumeDofs() = 0;
  
  /** Return number of d.o.f. per face. */
  virtual int getNfaceDofs()   = 0;
  
  /** Return number of d.o.f. per edge. */
  virtual int getNedgeDofs()   = 0;
  
  /** Return VTK type of volume element. */
  virtual int getVtk3Dtype() = 0;
  
  /** Return VTK type of face element. */
  virtual int getVtk2Dtype() = 0;
  
  /** Return VTK type of edge element. */
  virtual int getVtk1Dtype() = 0;
  
  /**
    Return global index of volume node in VTK ordering.
    @param volume Volume number.
    @param node   Local node number (in VTK ordering).
  */
  virtual int getVtkVolumeNode(int volume, int node) { return getVolumeNode(volume, node); }
};


/** Linear mesh. */
class LinMesh : public Mesh
{
 private:
  
  int addEdge(int n1, int n2, int nphys = 0, int ngeom = 0, int npart = 0);
  int addFace(int n1, int n2, int n3, int nphys = 0, int ngeom = 0, int npart = 0);
  
 public:

  LinMesh() { type = FE_LIN; };

  void writeGMSH(string fname) {}
  void findChildrenParents(bool create_new);
  
  int getNvolumeDofs() { return 4; }
  int getNfaceDofs()   { return 3; }
  int getNedgeDofs()   { return 2; }
  
  int getVtk3Dtype() { return VTK_TETRA; }
  int getVtk2Dtype() { return VTK_TRIANGLE; }
  int getVtk1Dtype() { return VTK_LINE; }
};


class QuadMesh : public Mesh
{
 private:
  int addEdge(int n1, int n2, int n3, int nphys = 0, int ngeom = 0, int npart = 0);
  
 public:

  QuadMesh() { type = FE_QUAD; };
 
  void writeGMSH(string fname);
  void findChildrenParents(bool create_new);
  
  int getNvolumeDofs() { return 10; }
  int getNfaceDofs()   { return 6; }
  int getNedgeDofs()   { return 3; }
  
  int getVtk3Dtype() { return VTK_QUAD_TETRA; }
  int getVtk2Dtype() { return VTK_QUAD_TRIANGLE; }
  int getVtk1Dtype() { return VTK_QUAD_EDGE; }
  int getVtkVolumeNode(int volume, int node) { return getVolumeNode(volume, VTK_QUAD_TETRA_NODE_ORDER[node]); }
};





#endif

